!**************************************************************************************************************************************
! This file contains stuff to calculate the linear power spectrum.
! You can choose from several transfer functions: BBKS, B & E, E & H (see below),
! and sigma_8 or CMB normalization.
!
! Declare an lps type:
!
! type(tLPS) :: lps
!
! and set it up with new:
!
! call new(lps, [universe], [mass_nu], [spectralIndex], [run], [transferFunction], [normalization])
!
! Optional arguments: [universe] is a type(tUniverse) that contains the other cosmological parameters not associated
!                                with the linear power spectrum. See base/universe.F90. If you don't supply one, it will use
!                                the default parameters.
!
!                     [spectralIndex] to change the value of n_s from the default in base/parameters.F90
!
!                     [run] change alpha from default in base/parameters.F90
!
!                     [transferFunction] pick your transfer function using a two-character string:
!                                        'BB' --- BBKS (default) good choice unless you want massive neutrinos.
!                                        'BE' --- Bond & Efstathiou.
!                                        'EH' --- Eisenstein & Hu's transfer function including massive neutrinos.
!                                        'BA' --- Eisenstein & Hu's transfer function including baryon oscillations.
!                                        'NW' --- Eisenstein & Hu's 'no-wiggles' transfer function, which is a smoothed version of 'BA'
!
!                     [normalization] 'S8' --- sigma_8 normalization (sigma_8 value is in the universe type)
!                                     'CB' --- CMB normalization using the WMAP scalar amplitude parameter A_s (also in universe type)
!
! Calculate the linear power spectrum using the function P_L(lps, k, a)
!
! Also there is the function smoothedVariance(lps, R, a, window), for window = 'TH' (top hat) or 'GA' (Gaussian).
!
!**************************************************************************************************************************************
module mLPS

        use mPrecision
        use mObject
        use mUniverse

        implicit none

        real(D), parameter :: N_NU = 3.0_D, &              ! Parameters for
                              T_CMB = 2.726_D, &           ! Eisenstein & Hu
                              THETA_CMB = T_CMB / 2.7_D, & ! transfer function
                              BARYON_T_CMB = 2.728_D, &
                              BARYON_THETA_CMB = BARYON_T_CMB / 2.7_D, &
                              K_PIVOT = 0.05_D, & ! 1 / Mpc
                              OPTICAL_DEPTH = 0.17_D, &
                              Z_MATTER_DOMINATED = 1088.0_D, &
                              A_MATTER_DOMINATED = 1.0_D / (1.0_D + Z_MATTER_DOMINATED)

        interface new
                module procedure new_lps
        end interface new

        interface write
                module procedure write_lps
        end interface write

        interface getUniverse
                module procedure getUniverse_lps
        end interface getUniverse

        interface setUniverse
                module procedure setUniverse_lps
        end interface setUniverse

        interface getH
                module procedure getH_lps
        end interface getH

        interface getOmega_m0
                module procedure getOmega_m0_lps
        end interface getOmega_m0

        interface getOmega_v0
                module procedure getOmega_v0_lps
        end interface getOmega_v0

        interface getOmega_b0
                module procedure getOmega_b0_lps
        end interface getOmega_b0

        interface getOmega_b
                module procedure getOmega_b_lps
        end interface getOmega_b

        interface setOmega_b
                module procedure setOmega_b_lps
        end interface setOmega_b

        interface getOmega_d
                module procedure getOmega_d_lps
        end interface getOmega_d

        interface setOmega_d
                module procedure setOmega_d_lps
        end interface setOmega_d

        interface getOmega_de
                module procedure getOmega_de_lps
        end interface getOmega_de

        interface setOmega_de
                module procedure setOmega_de_lps
        end interface setOmega_de

        interface omegas
                module procedure omegas_lps
        end interface omegas

        interface getW_0
                module procedure getW_0_lps
        end interface getW_0

        interface setW_0
                module procedure setW_0_lps
        end interface setW_0

        interface getW_a
                module procedure getW_a_lps
        end interface getW_a

        interface setW_a
                module procedure setW_a_lps
        end interface setW_a

        interface setB_1
                module procedure setB_1_lps
        end interface setB_1

        interface setB_2
                module procedure setB_2_lps
        end interface setB_2

        interface getB_1
                module procedure getB_1_lps
        end interface getB_1

        interface getB_2
                module procedure getB_2_lps
        end interface getB_2

        interface getShapeParameter
                module procedure getShapeParameter_lps
        end interface getShapeParameter

        interface P_L
                module procedure P_L_lps
        end interface P_L

        interface smoothedVariance
                module procedure smoothedVariance_lps
        end interface smoothedVariance

        interface getSpectralIndex
                module procedure getSpectralIndex_lps
        end interface getSpectralIndex

        interface setSpectralIndex
                module procedure setSpectralIndex_lps
        end interface setSpectralIndex

        interface getRun
                module procedure getRun_lps
        end interface getRun

        interface setRun
                module procedure setRun_lps
        end interface setRun

        interface getMass_nu
                module procedure getMass_nu_lps
        end interface getMass_nu

        interface getA_s
                module procedure getA_s_lps
        end interface getA_s

        interface getSigma_8
                module procedure getSigma_8_lps
        end interface getSigma_8

        private
        public :: K_PIVOT
        public :: new, write
        public :: getUniverse, setUniverse
        public :: getH
        public :: getW_0, setW_0
        public :: getW_a, setW_a
        public :: getA_s
        public :: getB_1, setB_1
        public :: getB_2, setB_2
!REMOVE        public :: norm, delta_zeta
        public :: norm
        public :: omegas
        public :: getSpectralIndex, setSpectralIndex
        public :: getRun, setRun
        public :: getMass_nu
        public :: transferFunction, baryonTransferFunction, TF_nowiggles
        public :: P_L
        public :: smoothedVariance, variance, k_max
        public :: topHatWindow

contains

        subroutine new_lps(lps, universe, mass_nu, spectralIndex, run, &
                                transferFunction, normalization)

                use mParameters

                type(tLPS), intent(out), optional :: lps
                type(tUniverse), intent(in), optional :: universe
                real(D), intent(in), optional :: mass_nu, &
                                                 spectralIndex, &
                                                 run
                character(len = TRANSFER_FUNCTION_LEN), parameter :: DEFAULT_TRANSFER_FUNCTION = 'BB'

                character(len = NORMALIZATION_LEN), parameter :: DEFAULT_NORMALIZATION = 'S8'

                character(len = TRANSFER_FUNCTION_LEN), intent(in), optional :: transferFunction

                character(len = NORMALIZATION_LEN), intent(in), optional :: normalization

                if( present(universe) ) then
                        lps%universe = universe
                else
                        call new(lps%universe)
                end if

                if( present(mass_nu) ) then
                        lps%mass_nu = mass_nu
                else
                        lps%mass_nu = DEFAULT_MASS_NU
                end if

                if( present(spectralIndex) ) then
                        lps%spectralIndex = spectralIndex
                else
                        lps%spectralIndex = DEFAULT_SPECTRAL_INDEX
                end if

                if( present(run) ) then
                        lps%run = run
                else
                        lps%run = DEFAULT_RUN
                end if

                if( present(transferFunction) ) then
                        lps%transferFunction = transferFunction
                else
                        lps%transferFunction = DEFAULT_TRANSFER_FUNCTION
                end if

                if( present(normalization) ) then
                        lps%normalization = normalization
                else
                        lps%normalization = DEFAULT_NORMALIZATION
                end if

                call normalize(lps)

        end subroutine new_lps

        subroutine write_lps(lps)

                type(tLPS), intent(in) :: lps

1               format(a, a)
2               format(a, f5.2)
3               format(a, f5.2, a)

                call write( lps%universe )

                write(6, 1) '# Transfer Function: ', lps%transferFunction
                write(6, 2) '# Spectral Index:', lps%spectralIndex
                write(6, 2) '# Run:', lps%run
                write(6, 3) '# mass_nu:', lps%mass_nu, ' eV'

        end subroutine write_lps

        type(tUniverse) pure function getUniverse_lps(lps)

                type(tLPS), intent(in) :: lps

                getUniverse_lps = lps%universe

        end function getUniverse_lps

        subroutine setUniverse_lps(lps, u)

                type(tLPS), intent(inout) :: lps
                type(tUniverse), intent(in) :: u

                lps%universe = u

        end subroutine setUniverse_lps

        real(D) function getH_lps(lps)

                type(tLPS), intent(in) :: lps

                getH_lps = getH( getUniverse(lps) )

        end function getH_lps

        real(D) function getOmega_m0_lps(lps)

                type(tLPS), intent(in) :: lps

                getOmega_m0_lps = getOmega_m0( getUniverse(lps) )

        end function getOmega_m0_lps

        real(D) function getOmega_v0_lps(lps)

                type(tLPS), intent(in) :: lps

                getOmega_v0_lps = getOmega_v0( getUniverse(lps) )

        end function getOmega_v0_lps

        real(D) function getOmega_b0_lps(lps)

                type(tLPS), intent(in) :: lps

                getOmega_b0_lps = getOmega_b0( getUniverse(lps) )

        end function getOmega_b0_lps

        real(D) function getOmega_b_lps(lps)

                type(tLPS), intent(in) :: lps

                getOmega_b_lps = getOmega_b( getUniverse(lps) )

        end function getOmega_b_lps

        subroutine setOmega_b_lps(lps, omega_b)

                type(tLPS), intent(inout) :: lps
                real(D), intent(in) :: omega_b

                type(tUniverse) :: u

                u = getUniverse(lps)
                call setOmega_b(u, omega_b)
                call setUniverse(lps, u)
                call normalize(lps)

        end subroutine setOmega_b_lps

        real(D) function getOmega_d_lps(lps)

                type(tLPS), intent(in) :: lps

                getOmega_d_lps = getOmega_d( getUniverse(lps) )

        end function getOmega_d_lps

        subroutine setOmega_d_lps(lps, omega_d)

                type(tLPS), intent(inout) :: lps
                real(D), intent(in) :: omega_d

                type(tUniverse) :: u

                u = getUniverse(lps)
                call setOmega_d(u, omega_d)
                call setUniverse(lps, u)
                call normalize(lps)

        end subroutine setOmega_d_lps

        real(D) function getOmega_de_lps(lps)

                type(tLPS), intent(in) :: lps

                getOmega_de_lps = getOmega_de( getUniverse(lps) )

        end function getOmega_de_lps

        subroutine setOmega_de_lps(lps, omega_de)

                type(tLPS), intent(inout) :: lps
                real(D), intent(in) :: omega_de

                type(tUniverse) :: u

                u = getUniverse(lps)
                call setOmega_de(u, omega_de)
                call setUniverse(lps, u)
                call normalize(lps)

        end subroutine setOmega_de_lps

        subroutine omegas_lps(lps, a, omega_m, omega_v)

                type(tLPS), intent(in) :: lps
                real(D) :: a, omega_m, omega_v

                call omegas( getUniverse(lps), a, omega_m, omega_v)

        end subroutine omegas_lps

        real(D) function getW_0_lps(lps)

                type(tLPS), intent(in) :: lps

                getW_0_lps = getW_0( getUniverse(lps) )

        end function getW_0_lps

        subroutine setW_0_lps(lps, w_0)

                type(tLPS), intent(inout) :: lps
                real(D), intent(in) :: w_0

                type(tUniverse) :: u

                u = getUniverse(lps)
                call setW_0(u, w_0)
                call setUniverse(lps, u)
                call normalize(lps)

        end subroutine setW_0_lps

        real(D) function getW_a_lps(lps)

                type(tLPS), intent(in) :: lps

                getW_a_lps = getW_a( getUniverse(lps) )

        end function getW_a_lps

        subroutine setW_a_lps(lps, w_a)

                type(tLPS), intent(inout) :: lps
                real(D), intent(in) :: w_a

                type(tUniverse) :: u

                u = getUniverse(lps)
                call setW_a(u, w_a)
                call setUniverse(lps, u)
                call normalize(lps)

        end subroutine setW_a_lps

        real(D) function getB_1_lps(lps)

                type(tLPS), intent(in) :: lps

                getB_1_lps = getB_1( getUniverse(lps) )

        end function getB_1_lps

        subroutine setB_1_lps(lps, b_1)

                type(tLPS), intent(inout) :: lps
                real(D), intent(in) :: b_1

                type(tUniverse) :: u

                u = getUniverse(lps)
                call setB_1(u, b_1)
                call setUniverse(lps, u)

        end subroutine setB_1_lps

        real(D) function getB_2_lps(lps)
                type(tLPS), intent(in) :: lps
                getB_2_lps = getB_2( getUniverse(lps) )
        end function getB_2_lps

        subroutine setB_2_lps(lps, b_2)

                type(tLPS), intent(inout) :: lps
                real(D), intent(in) :: b_2

                type(tUniverse) :: u

                u = getUniverse(lps)
                call setB_2(u, b_2)
                call setUniverse(lps, u)

        end subroutine setB_2_lps

        real(D) pure function getSpectralIndex_lps(lps)
                type(tLPS), intent(in) :: lps
                getSpectralIndex_lps = lps%spectralIndex
        end function getSpectralIndex_lps

        subroutine setSpectralIndex_lps(lps, index)
                type(tLPS), intent(inout) :: lps
                real(d), intent(in) :: index
                lps%spectralIndex = index
                call normalize(lps)
        end subroutine setSpectralIndex_lps

        real(D) pure function getRun_lps(lps)
                type(tLPS), intent(in) :: lps
                getRun_lps = lps%run
        end function getRun_lps

        subroutine setRun_lps(lps, run)
                type(tLPS), intent(inout) :: lps
                real(D), intent(in) :: run
                lps%run = run
                call normalize(lps)
        end subroutine setRun_lps

        real(D) pure function getMass_nu_lps(lps)
                type(tLPS), intent(in) :: lps
                getMass_nu_lps = lps%mass_nu
        end function getMass_nu_lps

        real(D) pure function getA_s_lps(lps)
                type(tLPS), intent(in) :: lps
                getA_s_lps = getA_s( getUniverse(lps) )
        end function getA_s_lps

        real(D) function getSigma_8_lps(lps)
                type(tLPS), intent(in) :: lps
                getSigma_8_lps = getSigma_8( getUniverse(lps) )
        end function getSigma_8_lps

        real(D) pure function getShapeParameter_lps(lps)
                type(tLPS), intent(in) :: lps
                getShapeParameter_lps = getShapeParameter( getUniverse(lps) )
        end function getShapeParameter_lps

        character(len = TRANSFER_FUNCTION_LEN) pure function transferFunction(lps)
                type(tLPS), intent(in) :: lps
                transferFunction = lps%transferFunction
        end function transferFunction
        
        subroutine normalize(lps)

                use mConstants

                type(tLPS), intent(inout) :: lps

                select case(lps%normalization)
                case default
                        write(0, *) 'ERROR: Improper normalization selection.'
                        stop
                case('S8') ! sigma_8 normalization
                        lps%norm = 1.0_D
                        lps%norm = getSigma_8(lps) ** 2 / smoothedVariance( lps, 8.0_D / getH(lps), 1.0_D, 'TH')
                case('CB') ! COBE normalization
                        lps%norm = (4.0_D / 25.0_D) &
!REMOVE                                   * delta_zeta() ** 2 &
                                   * getA_s(lps) / 1.84e4_D ** 2 &
                                   * (c * K_PIVOT &
                                      / ( 1e5_D * getH(lps) )) ** 4 &
                                   * ( A_MATTER_DOMINATED &
                                       / growthFactor(lps%universe, &
                                                    A_MATTER_DOMINATED) ) ** 2 &
                                   / getOmega_m0(lps) ** 2 &
                                   * (2.0_D * PI ** 2 / K_PIVOT ** 3) &
                                   / K_PIVOT ** getSpectralIndex(lps)

                end select
        end subroutine normalize

        real(D) function norm(lps)
                type(tLPS), intent(in) :: lps

                norm = lps%norm
        end function norm

!        real(D) function delta_zeta()
!                delta_zeta = 5.07e-5_D * exp( -(0.17_D - OPTICAL_DEPTH) )
!        end function delta_zeta

        real(D) function P_L_lps(lps, k, a)
                use mMisc
                type(tLPS), intent(in) :: lps
                real(D), intent(in) :: k, a ! wavenumber in units of 1 / Mpc

                ! Safety net:
!                call assert(lps%transferFunction == 'BA', 'wrong transfer function')

                select case (lps%transferFunction)
                case default
                        write(0, *) 'ERROR: I don''t know which linear power'
                        write(0, *) '       spectrum you want.'
                        write(0, *) lps%transferFunction
                        stop
                case('BE')
                        P_L_lps = BE_linearPowerSpectrum(lps%universe, k)
                case('BB')
                        P_L_lps = TFLinearPowerSpectrum(lps, k, a)
                case('EH')
                        P_L_lps = TFLinearPowerSpectrum(lps, k, a)
                case('BA')
                        P_L_lps = TFLinearPowerSpectrum(lps, k, a)
                case('NW')
                        P_L_lps = TFLinearPowerSpectrum(lps, k, a)
                end select
        end function P_L_lps

        !**********************************************************
        ! Calculates the linear power spectrum as P_L(k)=k^n*T^2(k)
        ! using the transfer function specified in parameters.in
        !**********************************************************
        real(D) function TFLinearPowerSpectrum(lps, k, a)
                type(tLPS), intent(in) :: lps
                real(D), intent(in) :: k, a ! 1 / Mpc
                real(D) :: transferFunction, kOverKPivot
                select case(lps%transferFunction)
                case default
                        write(0, *) 'ERROR: You are executing some code that'
                        write(0, *) '       should be impossible to reach.'
                        write(0, *) '       I conclude, therefore, that you do'
                        write(0, *) '       not exist!'
                        stop
                case('BB')
                        transferFunction = BBKSTransferFunction(lps, k)
                case('EH')
                        transferFunction = EHTransferFunction(lps, k, a)
                case('BA')
                        transferFunction = baryonTransferFunction(lps, k)
                case('NW')
                        transferFunction = TF_nowiggles(lps, k)
                end select
                if(lps%run == 0.0_D) then
                        TFLinearPowerSpectrum = lps%norm &
                                * k ** lps%spectralIndex &
                                * transferFunction * transferFunction
                else
                        kOverKPivot = k / K_PIVOT
                        TFLinearPowerSpectrum = lps%norm &
                            * k ** lps%spectralIndex &
                            * kOverKPivot ** ( 0.5_D * lps%run &
                                               * log(kOverKPivot) ) &
                            * transferFunction * transferFunction
                end if
        end function TFLinearPowerSpectrum

        !*****************************************************************
        ! Calculates the transfer function T(k) of Bardeen et. al. (BBKS).
        !*****************************************************************
        real(D) function BBKSTransferFunction(lps, k)

                type(tLPS), intent(in) :: lps
                real(D), intent(in) :: k
                real(D), parameter :: C1 = 3.89_D, &
                                      C2 = 16.1_D ** 2, &
                                      C3 = 5.46_D ** 3, &
                                      C4 = 6.71_D ** 4
                real(D) :: q, factor
                q = k / ( getShapeParameter(lps) * getH(lps) &
!                    * exp( -omega_b0(lps%universe) &
                    * exp( -getOmega_b0(lps) &
                    * (1.0_D + sqrt( 2.0_D * getH(lps) ) &
                    / getOmega_m0(lps)) ) )
                factor = 2.34_D * q
                BBKSTransferFunction = log(1.0_D + factor) / factor &
                    * ( 1.0_D + q * (C1 + q * ( C2 + q * (C3 + q * C4) )) ) &
                        ** (-0.25_D)

        end function BBKSTransferFunction

        real(D) pure function BE_LinearPowerSpectrum(u, k)

                use mConstants

                type(tUniverse), intent(in) :: u
                real(D), intent(in) :: k ! wavenumber in 1 / Mpc
                real(D) :: keff, q, q8, tk, tk8, Gamma
                Gamma = getShapeParameter(u)
                keff = log(Gamma / 0.36_D)
                keff = 0.172_D + 0.011_D * keff * keff
                q = 1e-20_D + k / ( Gamma * getH(u) )
                q8 = 1e-20_D + keff / Gamma
                tk = BEFit(q)
                tk8 = BEFit(q8)

                ! The formula above fits Delta(k):=k**3/(2PI**2) P(k)
                ! I strip off the k**3/(2PI**2) below and return P(k)
                BE_LinearPowerSpectrum = q / q8
                BE_LinearPowerSpectrum = PI * getSigma_8(u) * tk / tk8 &
                    * BE_LinearPowerSpectrum * BE_LinearPowerSpectrum
                BE_LinearPowerSpectrum = 2.0_D/ k ** 3 &
                    * BE_LinearPowerSpectrum * BE_LinearPowerSpectrum
        end function BE_LinearPowerSpectrum

        real(D) pure function BEFit(q)
                real(D), intent(in) :: q
                BEFit = 1.0_D / (1.0_D + &
                    ( 6.4_D * q + (3.0_D * q) ** 1.5_D &
                                + (1.7_D * q) ** 2) ** 1.13_D ) &
                    ** (1.0_D / 1.13_D)
        end function BEFit

        real(D) function smoothedVariance_lps(lps, R, a, window)
                use mIntegration
                type(tLPS), intent(in) :: lps
                real(D), intent(in)  :: R, a ! smoothing scale in Mpc
                character(len = INTEGRAL_CHAR_LEN), intent(in) :: window
                real(D), parameter :: LL = 1e-5_D, & ! Lower cutoff
                                      UL = 1e4_D     ! Upper cutoff
!***** I've estimated the error from introducing these cutoffs to be less than
!***** 0.0000001%, which you can check using the code commented out with !err.
#ifdef DEV
                real(D), parameter :: EPS = epsilon(1.0_R4)
#else ! DEV
!PRECISION       real(D), parameter :: EPS = epsilon(1.0_d)
                real(D), parameter :: EPS = epsilon(1.0_R4)
#endif ! DEV
!***************************************************************************
!err             real(d),parameter :: UL2=1e10_d ! You can use UL2 to see if
!err             claimed error from UL is really true. Adds integral from
!err             UL->UL2 and computes error from neglecting this part. Do
!err             something similar for LL.
!***************************************************************************
                type(tIntegral) :: svIntegral

                if(window /= 'TH' .and. window /= 'GA') then
                        write(0, *) 'ERROR: Window function not properly'
                        write(0, *) '       specified in SmoothedVariance.'
                        write(0, *) window
                        stop
                end if

                call new(svIntegral, LL, UL, EPS, changeVariables = 'EXP', &
                    params = (/ R, a /), char = window, lps = lps)
                call integrate(svIntegral, svIntegrand, smoothedVariance_lps)

!                if( .not. converged(svIntegral) ) then
!                        write(0, *) 'INFO: SmoothedVariance integration failed'
!                        write(0, *) '      to converge.'
!                end if
                call checkConvergence(svIntegral, 'SmoothedVariance integration failed to converge.')
!                call write(lps)
!                write(0, *) R, a, window

!err             call integrate(svIntegrand, UL, UL2, R, info, 'N', 'N')
!err             write(6, *) 'Relative error: ', UL / (SV + UL)
        end function smoothedVariance_lps

        !*****************************************************
        ! The integrand for calculating the smoothed variance.
        !*****************************************************
        real(D) function SVIntegrand(integral, k)
                use mConstants
                use mIntegration
                type(tIntegral), intent(in) :: integral
                real(D), intent(in) :: k ! wavenumber
                character(len = INTEGRAL_CHAR_LEN) :: windowSelector
                real(D) :: R, a, kR, window
                type(tLPS) :: lps
                call getParams(integral, R, a)
                windowSelector = getChar(integral)
                lps = getLPS(integral)
                !print *,'SVIntegrand: k = ',k

                kR=k*R
                !print *,'kR = ',kR
                select case(windowSelector)
                case default
                        write(0, *) 'ERROR: Window function not properly'
                        write(0, *) '       specified in svIntegrand.'
                        write(0, *) windowSelector
                        stop

                case('TH')
                        window = topHatWindow(kR)
                        window = k * window
                case('GA')
                        window = gaussianWindow(kR)
                        window = k * window
                end select
                !print *,'window = ',window,' PL = ',P_L(lps,k,a)
                SVIntegrand = window * window / (2.0_D * PI * PI) &
                    * P_L(lps, k, a)
                !print *,'SVIntegrand(',k,') = ',SVIntegrand
        end function svIntegrand

        real(D) function variance(lps, k_max, a)

                use mConstants
                use mIntegration

                type(tLPS), intent(in) :: lps
                real(D), intent(in) :: k_max, a
                real(D), parameter :: LL = 1e-5_D
                real(D), parameter :: EPS = 1.0e-4_D ! DEREK
                integer :: i
                type(tIntegral) :: integral
                call new(integral, LL, k_max, eps = EPS, &
                         changeVariables = 'EXP', lps = lps, params = (/ a /) )
                call integrate(integral, varianceIntegrand, variance)
                call checkConvergence(integral, 'INFO: variance integral failed to converge')
                variance = growthFactor( getUniverse( getLPS(integral) ), a) ** 2 * variance / (2.0_D * PI ** 2)

        end function variance

        real(D) function varianceIntegrand(integral, k)

                use mIntegration

                type(tIntegral), intent(in) :: integral
                real(D), intent(in) :: k
                real(D) :: a
                a = getParam(integral)
                varianceIntegrand = k ** 2 * P_L( getLPS(integral), k, a )

        end function varianceIntegrand

        !*****************************************************************
        ! Returns wavenumber where variance goes non-linear, no smoothing.
        !*****************************************************************
        real(D) function k_max(lps, a)

                use mRoot

                type(tLPS), intent(in) :: lps
                real(D), intent(in) :: a
                real(D), parameter :: ROOT_MIN = 1e-2_D, &
                                      ROOT_MAX = 1e2_D

                type(tRoot) :: root

                call new(root, ROOT_MIN, ROOT_MAX, params = (/ a /), lps = lps)
!                call findRoot(root, k_max_root, k_max)
                if( .not.findRoot(root, k_max_root, k_max) ) stop

        end function k_max

        real(D) function k_max_root(root, k)

                use mRoot

                type(tRoot), intent(inout) :: root
                real(D), intent(in) :: k

                real(D) :: a

                a = getParam(root)
                k_max_root = variance( getLPS(root), k, a ) - 1.0_D

        end function k_max_root

        real(D) pure function topHatWindow(kR)
                real(D), intent(in) :: kR
                topHatWindow = 3.0_D * ( sin(kR) - kR * cos(kR) ) / kR ** 3
        end function topHatWindow

        real(D) pure function gaussianWindow(kR)
                real(D), intent(in) :: kR
                gaussianWindow = exp(-0.5_D * kR * kR)
        end function gaussianWindow

        !*****************************************************************
        ! Eisenstein and Hu transfer function including neutrinos. This is
        ! code from Wayne Hu's webpage, modified to fit in with my code.
        !*****************************************************************
        real(D) function EHTransferFunction(lps, k, a)

                type(tLPS), intent(in) :: lps
                real(D), intent(in) :: k, a  ! 1/Mpc
                real(D) :: f_nu, f_baryon, &
                    omhh, kmax, kmin, dd_cb, dd_cbnu, dd0, T_master, obhh, &
                    z_equality, y_d, k_equality, r_drag, sound_horizon, &
                    alpha_nu, beta_c
                integer :: numk, i
                omhh = getOmega_m0(lps) * getH(lps) * getH(lps)
                ! f_nu from Eisenstein & Hu 1999:
                f_nu = lps%mass_nu * N_NU / (91.5_D * omhh)
                f_baryon = getOmega_b0(lps) / getOmega_m0(lps)
                if (f_nu .eq. 0.0_D) f_nu = 1e-10_D
                if (f_baryon .eq. 0.0_D) f_baryon = 1e-10_D

                T_CMB_sanity: if(T_CMB <= 0.0_D) then
                        write(0, *) 'Negative temperature requested!'
                        stop
                end if T_CMB_sanity
                N_NU_sanity: if (N_NU < 0.0_D) then
                        write(0, *) 'I can''t use a negative N_NU!'
                        stop
                end IF N_NU_sanity


                call TFset_parameters(omhh, f_baryon, f_nu, &
                    z_equality, y_d, sound_horizon, alpha_nu, beta_c)

                call growth(lps, k, a, omhh, f_nu, z_equality, &
                            DD_cb, DD_cbnu, DD0)

                T_master = TF_master(k, omhh, f_nu, alpha_nu, beta_c, &
                    sound_horizon)

                EHTransferFunction = T_master * DD_cbnu

        end function EHTransferFunction

        pure subroutine TFset_parameters(omhh,f_baryon,f_nu &
            ,z_equality,y_d,sound_horizon,alpha_nu,beta_c)

                real(D), intent(in) :: omhh, f_baryon, f_nu
                real(D), intent(out) :: z_equality, y_d, sound_horizon, &
                                        alpha_nu, beta_c
                real(D) :: obhh, k_equality, z_drag, &
                           r_drag, r_equality, P_c, P_cb, f_c, &
                           f_nub, f_cb

                obhh = omhh * f_baryon

                z_equality = 2.50e4_D * omhh * THETA_CMB ** (-4) - 1.0_D
                k_equality = 0.0746_D * omhh * THETA_CMB ** (-2)

                z_drag = 0.313_D * omhh ** (-0.419_D) &
                    * ( 1.0_D + 0.607_D * omhh ** (0.674_D) )
                z_drag = 1.0_D + z_drag * obhh &
                    ** ( 0.238_D * omhh ** (0.223_D) )
                z_drag = 1291.0_D * omhh ** (0.251_D) &
                    / ( 1.0_D + 0.659_D * omhh ** (0.828_D) ) * z_drag

                y_d = (1.0_D + z_equality) / (1.0_D + z_drag)

                R_drag = 31.5_D * obhh * THETA_CMB ** (-4) * 1000.0_D &
                    / (1.0_D + z_drag)
                R_equality = 31.5_D * obhh * THETA_CMB ** (-4) &
                    * 1000.0_D / (1.0_D + z_equality)

                sound_horizon = 2.0_D / (3.0_D * k_equality) &
                    * sqrt(6.0_D / R_equality) &
                    * log(( sqrt(1.0_D + R_drag) + sqrt(R_drag + R_equality) ) &
                    / ( 1.0_D + sqrt(R_equality) ))

                p_c  = -(5.0_D - sqrt( 1.0_D + 24.0_D &
                    * (1.0_D - f_nu - f_baryon) )) / 4.0_D
                p_cb = -(5.0_D - sqrt( 1.0_D + 24.0_D * (1.0_D - f_nu) )) &
                    / 4.0_D
                f_c  = 1.0_D - f_nu - f_baryon
                f_cb = 1.0_D - f_nu
                f_nub =f_nu + f_baryon

                alpha_nu = (f_c / f_cb) * ( 2.0_D * (p_c + p_cb) + 5.0_D ) &
                    / (4.0_D * p_cb + 5.0_D)
                alpha_nu = alpha_nu &
                    * (1.0_D - 0.553_D * f_nub + 0.126_D * f_nub ** 3)
                alpha_nu = alpha_nu &
                    / (1.0_D - 0.193_D * sqrt(f_nu) + 0.169_D *f_nu)
                alpha_nu = alpha_nu * (1.0_D + y_d) ** (p_c - p_cb)
                alpha_nu = alpha_nu &
                    * ( 1.0_D + (p_cb - p_c) / 2.0_D &
                        * (1.0_D + 1.0_D / ( (4.0_D * p_c + 3.0_D) &
                        * (4.0_D * p_cb + 7.0_D) )) &
                        / (1.0_D + y_d) )
                beta_c = 1.0_D / (1.0_D - 0.949_D * f_nub)

        end subroutine TFset_parameters

        real(D) pure function TF_master(k,omhh,f_nu,alpha_nu &
                                       ,beta_c,sound_horizon)

                real(D), intent(in) :: k
                real(D), intent(in) :: omhh, f_nu, alpha_nu, beta_c, &
                                       sound_horizon
                real(D) :: q, gamma_eff, q_eff, q_nu

                q = k * THETA_CMB ** 2 / omhh

                gamma_eff = ( sqrt(alpha_nu) + (1.0_D - sqrt(alpha_nu) ) &
                    / ( 1.0_D + (0.43_D * k * sound_horizon) ** 4) )

                q_eff = q / gamma_eff
                TF_master = log( exp(1.0_D) + 1.84_D * beta_c * sqrt(alpha_nu) &
                                 * q_eff )
                TF_master = TF_master / (TF_master + q_eff ** 2 * &
                    ( 14.4_D + 325.0_D / (1.0_D + 60.5_D * q_eff ** 1.11_D) ))

                q_nu = 3.92_D * q * sqrt(N_NU / f_nu)
                TF_master = TF_master &
                    * (1.0_D + ( 1.2_D * f_nu ** (0.64_D) &
                        * N_NU ** (0.3_D + 0.6_D * f_nu) ) &
                   / ( q_nu ** (-1.6_D) + q_nu ** (0.8_D) ))

        end function TF_master

        subroutine growth(lps, k, a, omhh, f_nu, z_equality, &
                               DD_cb, DD_cbnu, DD0)

                type(tLPS), intent(in) :: lps
                real(D), intent(in) :: k, a, omhh, f_nu, z_equality
                real(D), intent(out) :: DD_cb, DD_cbnu, DD0
                real(D) :: q, y_fs, oz, olz, P_cb, bigD, om, ov

                om = getOmega_m0(lps)
                ov = getOmega_v0(lps)

                q = k * THETA_CMB ** 2 / omhh

                y_fs = 17.2_D * f_nu &
                    * ( 1.0_D + 0.488_D * f_nu ** (-7.0_D / 6.0_D) ) &
                    * (N_NU * q / f_nu) ** 2

                oz = om * a ** (-3) &
                    / ( ov + (1.0_D - ov - om) * a ** (-2) + om * a ** (-3) )
                olz = ov &
                    / ( ov + (1.0_D - ov - om) * a ** (-2) + om * a ** (-3) )

                bigD = (1.0_D + z_equality) * a * 5.0_D * oz / 2.0_D &
                    * ( oz ** (4.0_D / 7.0_D) - olz + (1.0_D + oz / 2.0_D) &
                    * (1.0_D + olz / 70.0_D) ) ** (-1)

                DD0= bigD / ( (1.0_D + z_equality) * 5.0_D * om / 2.0_D &
                    * ( om ** (4.0_D / 7.0_D) - ov + (1.0_D + om / 2.0_D) &
                    * (1.0_D + ov / 70.0_D) ) ** (-1) )

                p_cb = -(5.0_D - sqrt( 1.0_D + 24.0_D * (1.0_D - f_nu) )) &
                        / 4.0_D

                DD_cb = ( 1.0_D + (bigD / (1.0_D + y_fs)) ** (0.7_D) ) &
                          ** (-p_cb / 0.7_D) &
                        * bigD ** p_cb

                DD_cbnu = ( (1.0_D - f_nu) ** (-0.7_D / p_cb) &
                            + ( bigD / (1.0_D + y_fs) ) ** (0.7_D) ) &
                          ** (-p_cb / 0.7_D) * bigD ** p_cb

        end subroutine Growth

        !*****************************************************************
        ! Eisenstein and Hu transfer function including baryons. This is
        ! code from Wayne Hu's webpage, modified to fit in with my code.
	! Ref: Eisenstein & Hu, ApJ, 1998, 496, 605.
        !*****************************************************************
        real(D) function baryonTransferFunction(lps, k)

                use mConstants

                type(tLPS), intent(in) :: lps
                real(D), intent(in) :: k
                real(D) :: f_baryon, k_equality, sound_horizon, k_silk, &
                           alpha_c, beta_c, alpha_b, beta_b, beta_node, q, &
                           ks, tf_cdm, s_tilde, tf_baryon

                call baryonTFset_parameters(lps, f_baryon, k_equality, &
                    sound_horizon, k_silk, alpha_c, beta_c, alpha_b, beta_b, &
                    beta_node)

                q = k / (13.41_D * k_equality)
                ks = k * sound_horizon

!  Main Variables

                tf_cdm = 1.0_D / ( 1.0_D + (ks / 5.4_D) ** 4 )
                tf_cdm = tf_cdm * TF_pressureless(q, 1.0_D, beta_c) &
                         + (1.0_D - tf_cdm) &
                         * TF_pressureless(q, alpha_c, beta_c)

                s_tilde = sound_horizon &
                    / ( 1.0_D + (beta_node / ks) ** 3.0_D ) ** ONE_THIRD
                tf_baryon = TF_pressureless(q, 1.0_D, 1.0_D) &
                            / ( 1.0_D + (ks / 5.2_D) ** 2 )
                tf_baryon = tf_baryon &
                    + alpha_b / ( 1.0_D + (beta_b / ks) ** 3 ) &
                              * exp( -(k / k_silk) ** 1.4_D )
                tf_baryon = tf_baryon * ( sin(k * s_tilde) / (k * s_tilde) )

                baryonTransferFunction = f_baryon * tf_baryon &
                        + (1.0_D - f_baryon) * tf_cdm

        end function baryonTransferFunction

        subroutine baryonTFset_parameters(lps, f_baryon, k_equality, &
            sound_horizon, k_silk, alpha_c, beta_c, alpha_b, beta_b, beta_node)

                type(tLPS), intent(in) :: lps
                real(D), intent(out) :: f_baryon, k_equality, sound_horizon, &
                                        k_silk, alpha_c, beta_c, alpha_b, beta_b, &
                                        beta_node


                real(D) :: omhh, obhh, z_equality, z_drag, R_drag, r_equality, y

                ! Input omhh -- The density of CDM and baryons, in units of
                ! critical dens, multiplied by the square of the Hubble constant, in
                ! units of 100 km/s/Mpc
                !
                ! f_baryon -- The fraction of baryons to CDM */
                ! Tcmb -- The temperature of the CMB in Kelvin, 2.728(4) is COBE
                !         and is the default reached by inputing Tcmb=0 
                !         -- reset on output.
                
                ! Output nothing, but set many global variables in common block 
                ! GLOBALVARIABLES. You can access them yourself, if you want:
                !
                !	theta_cmb,	/* Tcmb in units of 2.7 K */ 
                !	z_equality,	/* Redshift of matter-radiation equality,
                !                      really 1+z */
                !	k_equality,	/* Scale of equality, in Mpc^-1 */
                !	z_drag,		/* Redshift of drag epoch */
                !	R_drag,		/* Photon-baryon ratio at drag epoch */
                !	R_equality,	/* Photon-baryon ratio at equality epoch */
                !	sound_horizon,	/* Sound horizon at drag epoch, in Mpc */
                !	k_silk,		/* Silk damping scale, in Mpc^-1 */
                !	alpha_c,	/* CDM suppression */
                !	beta_c,		/* CDM log shift */
                !	alpha_b,	/* Baryon suppression */
                !	beta_b,		/* Baryon envelope shift */
                
                omhh = getOmega_m0(lps) &
                    * getH(lps) * getH(lps)
                f_baryon = getOmega_b0(lps) / getOmega_m0(lps)

                ! Are inputs reasonable?
                if (f_baryon.le.0.0_D) f_baryon=1.0e-5_D

                ! Auxiliary variables
                obhh = omhh*f_baryon

                ! Main variables
                z_equality = 2.50e4_D * omhh * BARYON_THETA_CMB ** (-4) - 1.0_D
                k_equality = 0.0746_D * omhh * BARYON_THETA_CMB ** (-2)

                z_drag = 0.313_D * omhh ** (-0.419_D) &
                    * ( 1.0_D + 0.607_D * omhh ** (0.674_D) )
                z_drag = 1.0_D + z_drag * obhh ** ( 0.238_D * omhh ** (0.223_D) )
                z_drag = 1291.0_D * omhh ** (0.251_D) &
                    / (1.0_D + 0.659_D * omhh ** (0.828_D) ) * z_drag

                R_drag = 31.5_D * obhh * BARYON_THETA_CMB ** (-4) &
                    * 1000.0_D / (1.0_D + z_drag) 
                R_equality = 31.5_D * obhh * BARYON_THETA_CMB ** (-4) &
                    * 1000.0_D / (1.0_D + z_equality) 

                sound_horizon = 2.0_D / (3.0_D * k_equality) &
                    * sqrt(6.0_D / R_equality) &
                    * log(( sqrt(1.0_D + R_drag) &
                    + sqrt(R_drag + R_equality) ) &
                    / ( 1.0_D + sqrt(R_equality) ))

                k_silk = 1.6_D * obhh ** (0.52_D) * omhh ** (0.73_D) &
                    * ( 1.0_D + (10.4_D * omhh) ** (-0.95_D) )

                alpha_c = ((46.9_D * omhh) ** (0.670_D) &
                    * ( 1.0_D + (32.1_D * omhh) ** (-0.532_D) ))
                alpha_c = alpha_c ** (-f_baryon) 
                alpha_c = alpha_c * ((12.0_D * omhh) ** (0.424_D) &
                    * ( 1.0_D + (45.0_D * omhh) ** (-0.582_D) )&
                    ) ** (-f_baryon ** 3)


                beta_c = 0.944_D / ( 1.0_D + (458.0_D * omhh) ** (-0.708_D) )
                beta_c = 1.0_D + beta_c &
                    * ( (1.0_D - f_baryon) &
                    ** ( (0.395_D * omhh) &
                    ** (-0.0266_D) ) - 1.0_D )
                beta_c = 1.0_D / beta_c

                y = (1.0_D + z_equality) / (1.0_D + z_drag)
                alpha_b = y * ( -6.0_D * sqrt(1.0_D + y) + (2.0_D + 3.0_D * y) &
                    * log(( sqrt(1.0_D + y) + 1.0_D ) &
                    / ( sqrt(1.0_D + y) - 1.0_D )) )
                alpha_b = 2.07_D * k_equality * sound_horizon &
                    * (1.0_D + R_drag) ** (-0.75_D) * alpha_b


                beta_b = 0.5_D + f_baryon + (3.0_D - 2.0_D * f_baryon) &
                    * sqrt( (17.2_D * omhh) ** 2 + 1.0_D )

                beta_node = 8.41_D * omhh ** (0.435_D)

        end subroutine baryonTFset_parameters

        real(D) pure function TF_pressureless(q, a, b)
                use mConstants

                real(D), intent(in) :: q, a, b

                TF_pressureless = log(E + 1.8_D * b * q)
                TF_pressureless = TF_pressureless / (TF_pressureless + &
                    ( 14.2_D / a + 386.0_D / (1.0_D + 69.9_D * q ** 1.08_D) ) &
                    * q ** 2)

        end function TF_pressureless

        real(D) function TF_nowiggles(lps, k)

                type(tLPS), intent(in) :: lps
                real(D), intent(in) :: k
                real(D) :: omhh, f_baryon, q_eff, a

                omhh = getOmega_m0(lps) * getH(lps) * getH(lps)
                f_baryon = getOmega_b0(lps) / getOmega_m0(lps)

                ! Are inputs reasonable?
                if(f_baryon <= 0.0_D) f_baryon = 1.0e-5_D

                a = alpha_gamma(omhh, f_baryon)
                q_eff = k / omhh * (BARYON_THETA_CMB) ** 2
                q_eff = q_eff / (a + (1.0_D - a) &
                    / (1.0_D + (0.43_D * k &
                    * sound_horizon_fit(omhh, f_baryon)) ** 4))

                TF_nowiggles = TF_zerobaryon(q_eff)

        end function TF_nowiggles

        real(D) pure function alpha_gamma(omhh, f_baryon)

                real(D), intent(in) :: omhh, f_baryon

                alpha_gamma = 1.0_D - 0.328_D * log(431.0_D * omhh) * f_baryon &
                    + 0.38_D * log(22.3_D * omhh) * (f_baryon) ** 2

        end function alpha_gamma

        real(D) pure function sound_horizon_fit(omhh, f_baryon)

                real(D), intent(in) :: omhh, f_baryon
                real(D) :: obhh
                obhh = f_baryon * omhh
                sound_horizon_fit = 44.5_D * log(9.83_D / omhh) &
                    / sqrt( 1.0_D + 10.0_D * obhh ** (0.75_D) )

        end function sound_horizon_fit

        real(D) pure function TF_zerobaryon(q)

                use mConstants

                real(D),intent(in) :: q

                TF_zerobaryon = log(2.0_D * E + 1.8_D * q)
                TF_zerobaryon = TF_zerobaryon / ( TF_zerobaryon &
                    +(14.2_D + 731.0_D / (1.0_D + 62.5_D * q) ) * q ** 2)

        end function TF_zerobaryon

end module mLPS
